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ABSTRACT 


Pulse  Detonation  Engines  (PDE)  are  seen  to  be  the  next  generation  propulsion 
systems  due  to  enhaneed  thermodynamie  efficiencies  based  on  the  Humphrey  cycle.  One 
of  the  limitations  in  fielding  practical  designs  has  been  attributed  to  tube  diameters  not 
exceeding  5  inches  as  the  shock  wave  takes  a  long  run  distance  for  transition  to 
detonation,  thus  potentially  affecting  specific  thrust.  Novel  methods  via  imploding 
detonations  were  investigated  to  remove  such  limitations.  During  the  study,  a  practical 
computational  cell  size  was  first  determined  so  as  to  capture  the  required  physics  for 
transient  detonation  wave  propagation  using  a  Hydrogen-Air  test  case.  Through  a  grid 
sensitivity  analysis,  one-quarter  of  the  induction  length  was  found  sufficient  to  capture 
the  experimentally  observed  initial  wave  transients.  Test  case  models  utilizing  axi- 
symmetric  head-on  implosions  were  studied  in  order  to  understand  how  the  implosion 
process  reinforces  a  detonation  wave  as  it  expands.  This  in  effect  creates  localized 
overdriven  regions,  which  maintains  the  transition  process  to  full  detonation.  A 
parametric  study  was  also  performed  to  determine  the  extent  of  diameter  increase  and  it 
was  found  that  the  detonations  could  be  supported  with  no  change  in  run  distance  even 
when  the  tube  diameter  exceeds  5  inches. 
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I.  INTRODUCTION 


A.  MOTIVATION 

Pulse  Detonation  Engines  (PDEs)  have  tremendous  potential  to  replaee  eurrent 
state-of-the  art  propulsion  systems,  whieh  ean  operate  up  to  supersonie  eruise  speeds  for 
missiles,  fighter  jets,  unmanned  aerial  vehieles,  ete.  Much  of  this  potential  stems  from  the 
fact  that  PDEs  operates  on  the  detonation  mode  of  combustion,  which  has  a  higher 
thermodynamic  efficiency  than  current  propulsion  systems  based  upon  constant  pressure. 
However,  no  practical  design  has  been  fielded  to  date  and  one  key  limitation  is  due  to  the 
fact  that  detonation  can  only  be  supported  in  main  detonation  chambers  of  up  to  five 
inches  in  diameter  without  requiring  unreasonably  long  tube  lengths.  This  impairs  the 
specific  thrust  substantially  as  multiple  tubes  are  required  to  generate  thrust  required  by  a 
typical  cruise  missile.  As  a  result,  thrust  density  decreases  and  potential  cross  flow 
mixing  problems  from  the  effiux  of  five  separate  tubes  reduce  the  efficiency  even  further. 
To  increase  the  diameter  of  the  main  tube  and  yet  maintain  the  tube  length,  simulation 
studies  involving  novel  means  using  shock-focusing  techniques  via  toroidal  wave 
implosion  first  investigated  by  Jackson  et  al.  [1]  is  conducted.  This  reinforces  the  wave 
propagation  process  and  effectively  shortens  the  deflagration-to-detonation  run  distance 
such  that  a  detonation  wave  can  be  supported  for  larger  diameter  tubes. 

B,  BACKGROUND 

Pulse  Detonation  Engines  (PDEs)  operates  on  the  detonation  mode  of 
combustion,  which  can  be  best  approximated  by  the  Humphrey  (Constant-Volume)  cycle 
in  contrast  to  the  Brayton  (Constant-Pressure)  cycle  for  gas  turbine  engines,  as  illustrated 
in  Eigure  1.  The  amount  of  work  done  (Pressure- Volume)  by  the  Humphrey  cycle  as  heat 
energy  is  introduced  under  constant  volume  conditions  from  2  to  3  is  much  greater  as 
compared  to  the  Brayton  cycle  under  constant  pressure  conditions,  resulting  in  higher 
thermodynamic  efficiencies. 
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Figure  1,  Pressure-volume  plot  based  on  Humphrey  cycle  [From  2] 

In  order  to  create  detonation  shock  waves  that  can  propagate  at  supersonic  speeds, 
a  finite  run  length  distance  is  required  for  such  transitions  from  defiagration-to- 
detonation  known  as  a  DDT  process  for  a  typical  confined  tube  configuration.  As  tube 
diameters  become  larger,  a  longer  DDT  run  distance  is  required  in  order  to  support  the 
transition  to  a  fully  developed  detonation  wave.  This  limitation  results  in  multiple  tube 
configurations  for  a  practical  cruise  missile  design. 

A  novel  method  is  being  proposed  using  shock  focusing  techniques  via  toroidal 
wave  implosion  such  that  the  shock  wave  propagation  in  the  main  shock  tube  is 
reinforced  and  consequently,  the  resulting  pressures  generated  are  much  higher  than  the 
Chapman- Jouguet  (C-J)  conditions  defined  for  a  fully-developed  steady  state  wave.  The 
extreme  high  pressure  environment  favors  much  higher  reaction  rates  for  the  fuel- 
oxidizer  reactant  mixtures  and  this  will  allow  for  a  near  direct  generation  of  a  detonation 
(without  the  need  for  DDT)  as  the  fuel-air  mixtures  react  violently  in  order  to  support  the 
shock  front. 

In  order  to  explore  such  implosion  concepts  accurately,  it  is  imperative  to  capture 
the  full  physics  and  chemical  kinetics  behind  such  complex  and  highly  transient  fluid 
flows  governing  the  shock  wave  behavior  followed  by  chemical  reactions  based  on  a 
reduced  set  of  chemical  reaction  mechanisms.  Simulations  of  such  complex  flow  field 
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behavior  are  only  achievable  via  a  practical  computational  cell  mesh  topology  that  is 
sufficiently  small  enough  on  the  order  of  the  induction  length  to  describe  the  reaction 
zone  behind  the  propagating  shock  front.  Validation  of  such  simulation  studies  based  on 
well  known  reduced  kinetic  models  for  Hydrogen-Air  mixtures  can  be  conducted  via 
comparison  with  experimental  data  available  in  the  literature. 

The  required  computational  cell  size  for  such  high  fidelity  simulations  can  then 
be  further  extended  to  generic  hydrocarbon  fuel-air  mixtures  such  as  acetylene-air, 
ethylene-air,  propane-air  and  even  JP  10-air  mixtures  based  on  the  required  length  of  the 
induction  zone  supporting  the  shock  front.  Utilizing  the  practical  computational  cell  size, 
further  analysis  for  a  simulation  case  model  based  on  an  axi-symmetric  head-on  toroidal 
implosion  can  be  conducted  and  analyzed  for  similar  configuration  types. 

Parametric  studies  are  then  performed  to  investigate  the  extent  of  the  increase  in 
tube  diameters  for  a  fixed  practical  propulsion  tube  length  based  on  the  test  case  model. 

C.  OBJECTIVES 

The  research  objective  is  to  enhance  current  PDE  applications  by  carrying  out 
simulation  studies  to  investigate  the  extent  of  possible  combustor  diameter  increase  using 
shock  focusing  techniques  via  toroidal  implosions  with  modified  PDE  tube 
configurations  without  introducing  the  need  for  longer  tube  lengths. 

Simulation  studies  were  conducted  on  the  CED++  computational  platform  and 
carried  out  via  the  following  order: 

i)  Determination  of  a  practical  computational  cell  size  for  high  fidelity 
capturing  of  detonation  wave  propagation  in  a  simple  PDE  tube 
configuration  based  on  well  known  Hydrogen-Air  reduced  reaction 
models 

ii)  Simulation  and  analysis  of  head-on  axi-symmetric  toroidal  implosions  for 
Hydrogen-Air  mixture  based  on  the  required  computational  cell  size 
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iii)  Performance  of  parametric  studies  used  to  evaluate  the  extent  of  the 
diameter  increase  for  a  fixed  practical  tube  length  using  toroidal  implosion 
based  on  an  enhanced  PDE  design  configuration 
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II.  THEORY 


A.  PULSE  DETONATION  ENGINES 

Pulse  Detonation  Engines  are  essentially  intermittent  propulsion  systems  where 
thrust  is  provided  by  repetitive  triggering  of  detonation  waves  propagating  at  supersonie 
speeds  in  a  long  tube  enelosed  at  the  head  end.  The  proeess  takes  plaee  under  nearly 
eonstant  volume  eonditions  and  its  operating  eyele  ean  be  effeetively  described  by  three 
distinct  phases.  The  first  phase  consists  of  introducing  a  reactive  component  such  as 
hydrogen  fuel  and  an  oxidizing  component  such  as  air  into  the  main  combustion 
chamber.  Initiation  of  the  fuel-air  mixture  near  the  enclosed  portion  of  the  chamber 
occurs  in  the  second  phase  via  an  energy  sparking  source.  A  high  pressure  detonation 
wave  as  a  result  of  the  initiation  ensues  in  the  third  phase.  The  remaining  fuel-air  mixture 
in  the  tube  is  compressed  in  a  thin  reaction  zone  by  the  traversing  detonation  wave  and 
reacts  violently  under  extreme  high  pressure  conditions  to  feed  the  leading  edge  of  the 
shock  front  followed  by  an  expansion  zone  of  detonation  products.  The  last  phase 
consists  of  effectively  purging  the  gaseous  products  out  of  the  main  chamber  and  the 
entire  cycle  repeats.  Thrust  is  produced  as  the  gaseous  products  exit  out  of  the  main 
propulsion  tube.  Figures  2  to  4  describes  each  phase  of  a  PDF  cycle. 
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Figure  2.  PDE  cycle  -  fill  phase 
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Figure  4,  PDE  cycle  -  purge  phase 


B,  DETONATION  THEORY 

Detonation  or  deflagration  by  itself  is  an  exothermie  reaetion  eaused  by  transition 
of  unstable  reaetive  molecular  reactions  to  a  stable  state.  The  amount  of  energy  required 
to  trigger  such  reactions  is  known  as  the  Activation  Energy.  When  such  energy  levels  are 
reached,  the  original  reactant  molecules  come  apart  and  re-arrange  themselves  back  into  a 
more  stable  state  producing  shock  waves  as  a  result  of  a  sudden  sharp  increase  in 
thermodynamic  properties  due  to  the  exothermic  reaction. 

Combustion  of  fuel-air  mixtures  produces  either  a  deflagration  or  detonation 
wave.  Properties  such  as  pressure  (P),  temperature  (T)  and  density  (p)  of  the  reactant  and 
products  across  the  combustion  wave  front  change  discontinuously.  Figure  5  illustrates 
the  change  across  such  properties  for  a  constant  area  duct  where  the  flame  front  is 
envisioned  as  stationary. 
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Figure  5,  One-dimensional  combustion  wave  model  [From  3] 


In  order  to  relate  properties  aeross  a  one-dimensional  eombustion  wave  front,  the 
eorresponding  mass  eonservation,  momentum  eonservation,  energy  eonservation  and 
equation  of  state  must  be  invoked.  For  the  state  equation,  the  ideal  gas  law  serves  as  a 
good  approximation  to  model  the  properties  of  the  gases  during  eombustion.  These 
relations  assume  that  the  eombustion  wave  is  steady  state  and  the  proeess  aeross  the  wave 
front  is  adiabatie  and  remain  in  ehemieal  and  thermodynamie  equilibrium  [3]. 

Mass  Conservation: 


PlMj  P2^2 

Momentum  Conservation: 


P\  +  Awf  =  A  +  P2^l 


Energy  Conservation: 


2  2 
U,  Un 

CX+^  +  q  =  CP+^ 


(1) 

(2) 


(3) 


Equation  of  State  (Ideal  Gas  Law): 

P  =  pRT  (4) 

Speeifie  heat  relation: 


C 


p 


7 

y-\ 


R 


(5) 
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where: 


R  =  Gas  Constant 

q  =  Speeifie  heat  energy  added  to  the  system  via  the  eombustion  proeess 
Cp  =  Speeifie  heat  at  eonstant  pressure 
Y  =  Ratio  of  speeifie  heats 

In  order  to  gain  further  insight  into  all  possible  states  that  the  burned  gas  produets 
ean  exist,  a  Hugoniot  relation  must  be  determined.  A  Hugoniot  curve  is  basically  a  plot  of 

pressure  (Rj)  versus  specific  volume  (  — )  corresponding  to  the  burned  gas  products  for 

Pi 

a  fixed  heat  release  mass,  q  and  the  relation  can  be  derived  by  combining  equations  (2), 
(3)  and  (4)  to  obtain  the  given  Hugoniot  relation. 

Hugoniot  Relation: 


r  { Pi 

Y-Api  Pij 


1  1 

- 1 - 

P\  Pi 


(6) 


The  slope  of  the  Hugoniot  relation  indicates  the  velocity  of  the  detonation  wave.  The 
Rayleigh  Line  relation  is  obtained  by  combining  equations  (1)  and  (2). 

Rayleigh  Line  Relation: 


(V) 


The  Hugoniot  curve  is  illustrated  in  Figure  6.  From  the  curve,  the  origin  A 

represents  the  pressure  (F;)  and  specific  volume  (  — )  of  the  unbumed  reactants.  The 

Px 

two  tangent  points  intersecting  origin  A  are  defined  as  upper  and  lower  Chapman-Jouget 
(C-J)  points.  From  the  Hugoniot  curve,  the  C-J  points  can  be  obtained  via  the  slope  of  the 
Hugoniot  curve  as  follows: 
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C-J  relation; 


dPz  jP^-P,) 

' .  /  A  r  ^  1  ^ 


4  y  ] 

fill 

WP2) 

1 

1 

(8) 


The  C-J  relation  ean  also  be  obtained  via  differentiation  of  equation  (6)  with  respect  to 


dP. 


(P2-P,)- 


^  2/  ^ 
7-1, 


iv] 

^  2y  ^ 

1 

rj__  n 

WP2J 

I/-I7 

P2 

V.P2  Ay 

(9) 


The  speed  of  sound,  C2  for  the  burned  products  is  given  by  the  following  relation: 


C2 


(10) 


Equating  (8)  and  (9)  along  with  the  Rayleigh-Line  relation  in  (7)  gives  the  following: 


u 


2 

2 


(11) 


This  implies  that  the  particle  velocity  of  the  gases  leaving  the  detonation  wave  are 
traveling  at  sonic  velocity.  In  a  shock-frame  of  reference,  this  means  that  the  product 
gases  leaving  the  detonation  wave  are  leaving  at  the  sonic  velocity  of  the  products. 
Therefore  equilibrating  the  amount  of  energy  that  reaches  the  shock  front  from  the 
chemical  reactions  and  thus  allowing  a  sustained  steady-state  detonation  wave  to 
propagate  at  what  is  known  as  the  Chapman-Jouguet  velocity.  Should  this  not  be  the 
case,  detonation  would  undergo  an  unsteady  process  that  will  most  likely  asymptote  to 
the  C-J  condition. 

Five  distinct  regions  can  be  inhered  from  the  Hugioniot  curve,  however,  not  all 
regions  are  physically  valid.  Table  1  describes  the  validity  of  the  states  in  these  regions. 
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Figure  6,  Hugoniot  curve  for  fuel-air  mixtures  [From  3] 


Table  1,  Physical  validity  of  states  on  the  Hugoniot  curve 


Region 

Remarks 

I 

Region  I  is  defined  as  transient  states  where  the  combustion  waves  propagate 

as  strong  overdriven  detonations  above  the  upper  C-J  point.  Under  such 

conditions,  the  lead  shock  front  travels  at  a  much  faster  rate  than  the  product 

gases  resulting  in  a  larger  induction  zone.  The  induction  zone  will  eventually 

increase  to  such  an  extent  such  that  any  heat  release  generated  by  the 

combustion  of  the  fuel-oxidizer  reactants  will  not  affect  the  lead  shock  front, 

thus  slowing  it  down.  As  such,  strong  detonation  waves  will  always  decay 

back  to  the  upper  C-J  point. 

II 

Region  II  is  defined  as  transient  states  where  the  combustion  waves  propagate 
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Region 

Remarks 

as  weak  detonations  below  the  upper  C-J  point.  Under  sueh  eonditions,  the 

veloeity  of  the  gas  produets  travels  at  a  faster  rate  than  the  lead  shoek  front  and 

as  sueh  the  spatial  extent  of  the  induetion  zone  deereases.  The  same  heat 

release  energy  generated  from  the  eombustion  reaetion  starts  to  have  an 

inereasing  effeet  on  the  lead  shoek  front  due  to  the  deereasing  induction  zone 

length  and  the  additional  heat  release  supplied  starts  to  accelerate  the  lead 

shock  front  away  from  the  reaction  zone  back  towards  the  upper  C-J  point. 

III 

Region  III  represents  states  where  the  combustion  waves  propagate  as  weak 

deflagrations.  They  exist  as  expansion  or  rarefaction  waves  where  the  density 

of  the  burned  products  is  lower  than  that  corresponding  to  the  reactants.  They 

are  observed  experimentally  during  the  initial  deflagration-to-detonation 

process  at  ignition  and  the  gas  velocity  relative  to  the  wave  front  is  accelerated 

from  lower  to  higher  sub-sonic  speeds. 

IV 

Region  IV  represents  states  where  the  combustion  waves  propagate  as  strong 

deflagration  waves.  These  are  physically  impossible  states  where  the  gas 

velocity  relative  to  the  wave  front  must  be  accelerated  from  sub-sonic  to  super¬ 
sonic  speeds,  which  are  forbidden  for  a  wave  structure  propagating  in  a 

constant  area  duct. 

V 

The  solution  states  for  region  V  is  mathematically  impossible.  The  solution  to 

the  Rayleigh-Line  equation  (7)  results  in  an  imaginary  number  for  Ui. 

From  the  analysis  given  in  Table  1,  the  physieally  validity  of  regions  1  and  3  are 
most  likely  observed  during  a  detonation  proeess  and  this  has  been  eonfirmed  in  shoek 
tube  experiments  as  well  [3]. 

C.  DEFLAGRATION-TO-DETONATION  TRANSITION 

A  detonation  wave  requires  a  finite  run  distanee  for  the  shoek  wave  to  transit  from 


a  deflagration  to  a  detonation  wave  known  as  DDT.  During  initiation  of  a  fuel-oxidizer 
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mixture  in  a  PDE  tube,  a  deflagration  leading  shoek  wave  ensues  and  propagates  at  a 
wave  speed  eharaeteristie  of  the  reaetant  temperature  eonditions  followed  by  a  series  of 
eompression  waves  produeed  by  the  eombustion  proeess.  The  eompression  waves  further 
heats  up  the  region  behind  the  leading  shoek  wave  ereating  loealized  high  temperature 
regions  and  eonsequently  inerease  the  veloeity  of  the  eompression  waves.  The 
eompression  waves  eventually  eateh  up  with  the  leading  wave  front  eoaleseing  into  a 
detonation  wave  [3].  The  tube  eonfinement  as  well  as  the  turbulenee  generated  by  the 
motion  of  the  produet  gases  from  the  detonation  shook  front  leads  to  the  onset  of  “an 
explosion  in  an  explosion,”  whioh  produoes  an  overdriven  forward  shook  wave  into  the 
unbumed  reaotants  known  as  superdetonation  and  an  opposite  shook  wave  moving  into 
the  gas  produots  known  as  retonation  with  transverse  osoillating  shooks  in  between.  This 
finally  leads  to  a  self-sustained  steady  state  detonation  wave  propagating  at  C-J  veloeity. 
Figure  7  illustrates  this  phenomenon. 


Figure  7,  Transition  process  to  superdetonation  [From  3] 


The  DDT  run  distanoe  is  impraotioally  long  for  ourrent  applioations  and  therefore 
obstaoles  suoh  as  swept  rams  or  spirals  are  used  in  PDE  oonfigurations  to  shorten  the 
DDT  length  via  turbulent  flow  regimes  to  promote  higher  ehemieal  reaetion  kineties. 
However,  generation  of  high  pressure  losses  as  a  eonsequenee  of  suoh  obstaoles  degrades 
the  propulsive  effioienoy  tremendously. 
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D,  DETONATION  WAVE  STRUCTURE 


Sustainment  of  the  detonation  wave  ean  be  deseribed  simply  by  the  one¬ 
dimensional  model  developed  independently  by  Zeldovich,  von  Neumann  and  Boring 
[3].  This  is  known  as  the  ZND  model  illustrated  in  Figure  8.  The  model  consists  of  a 
shock  wave  moving  at  detonation  velocity  followed  by  a  rarefaction  region.  The  shock 
wave  region  is  extremely  thin  in  the  order  of  a  few  mean  free  paths  thick  of  molecules; 
therefore  limited  reactions  occur  in  this  region.  Across  the  shock  region,  there  is  a  sharp 
almost  discontinuous  increase  of  pressure,  temperature  and  density.  The  rarefaction 
region  is  comparatively  much  thicker  usually  on  the  order  of  1  cm  and  consists  of  an 
induction  and  reaction  zone.  The  thermodynamic  properties  in  the  induction  zone 
immediately  behind  the  shock  front  remain  almost  constant.  Reaction  rates  as  a  function 
of  temperature  begin  to  increase  slowly.  Once  the  induction  period  is  completed, 
temperatures  start  to  climb  sharply  with  extremely  high  heat  release  energies  caused  by 
the  violent  reactions  of  the  fuel-oxidizer  mixture  defined  by  the  reaction  zone,  which 
sustains  the  propagation  of  the  shock  wave.  The  gasdynamic  properties  reach  their 
equilibrium  states  at  the  end  of  the  reaction  zone.  Therefore,  the  induction  period  or 
induction  zone  is  extremely  important  in  order  to  capture  the  correct  thermodynamic 
properties  of  temperature,  pressure  and  density  for  an  accurate  simulation  of  a  detonation 
wave. 
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Figure  8,  Thermodynamic  properties  of  the  ZND  model  [From  3] 


E.  DETONATION  CELL  WIDTH 

Detonation  cell  widths  are  multi-dimensional  detonation  wave  effects 
experimentally  measured  by  soot  foil  imprint  techniques  where  a  detonation  wave  created 
by  a  typical  fuel-air  mixture  leaves  a  fish  scale  pattern  signature  on  soot  coated  aluminum 
sheet  lined  on  the  inner  surface  of  a  shock  tube.  Figure  9  illustrates  a  typical  soot  foil 
imprint  produced  by  a  detonating  shock  wave  from  JPlO-Air  mixture  where  JPIO  is  a 
liquid  hydrocarbon  fuel  (CioHie).  The  non-planar  detonation  shock  front  induced  by  the 
energy  release  arising  from  the  combustion  of  the  fuel-oxidizer  mixture  is  a  resultant  of 
the  Mach  stem,  incident  shock  and  reflected  shock,  which  is  a  necessary  condition  for  a 
self-sustained  detonation  front.  The  interaction  of  the  three  waves  produces  a  shear 
discontinuity  known  as  the  triple  point.  Figure  10  illustrates  such  a  detonation  front 
structure.  As  a  result,  the  waves  propagate  down  the  tube  leaving  the  characteristic  fish 
scale  pattern  signature.  The  transverse  distance  between  each  fish  scale  is  known  as  the 
cell  width,  which  is  regarded  as  a  fundamental  parameter  to  characterize  the  detonability 
of  a  gaseous  fuel-air  mixture.  A  more  reactive  fuel  such  as  hydrogen  as  compared  to  a 
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typical  hydrocarbon  fuel  such  as  ethylene  (C2H4)  will  eorrespondingly  have  a  smaller 
detonation  eell  width.  Similarly,  reaetion  of  hydrogen  in  the  presenee  of  a  less  reaetive 
oxidizer  sueh  as  air  with  respeet  to  pure  oxygen  will  have  larger  eell  widths. 


Figure  9  Soot  foil  imprint  produced  by  JPlO-Air  detonation  wave  [From  4] 


Figure  10,  Detonation  front  structure  with  triple  point  [From  4] 
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F. 


REDUCED  CHEMICAL  KINETIC  REACTION  MODELS 


Chemical  kinetics  governs  the  rate  at  which  the  species  concentration  of  a  fuel- 
oxidizer  mixture  is  transformed  into  its  product  species  and  vice  versa  depending  on 
pressure  and  temperature  reaction  conditions  as  well  as  complexity  of  the  involved 
species.  The  dependence  of  the  reaction  rate  constant,  A:  as  a  function  of  temperature,  T 
can  be  inferred  from  the  Arrhenius  Equation  [3]. 

Arrhenius  Equation: 


k  =  AT"  exp 


E 

RTy 


(12) 


The  pre-exponential  parameter  A  describes  the  frequency  of  collisions  for  the  species  that 
result  in  a  reaction  as  well  as  the  preferred  orientation  of  the  molecular  species  that  favor 
a  reaction.  The  exponent  n  expresses  the  temperature  dependence  of  the  pre-exponential 
parameter.  Activation  energy,  Ea  defines  the  minimum  amount  of  energy  required  for  a 
reaction  to  occur  and  R  is  the  universal  gas  constant. 

While  it  is  desired  to  include  a  detailed  chemical  reaction  model  governing  a 
reactive  fuel-oxidizer  mixture,  which  can  involve  up  to  hundreds  of  reactions  and  species 
under  all  operating  conditions,  it  is  computationally  demanding  especially  for  complex 
flows  involving  cellular  detonation  wave  structures  even  for  super-computers  of  today. 
To  ensure  that  computational  times  required  for  CED  simulations  to  be  tractable  under 
reasonable  time  frames,  it  is  imperative  to  reduce  detailed  chemical  mechanisms  into 
simpler  models  where  only  certain  chemical  reactions  are  significant  for  example  during 
a  detonation  process  and  yet  able  to  capture  the  physics  of  the  shock  wave  propagation 
with  correct  heat  release  rates  in  the  reaction  zone.  This  can  be  done  via  approximation  of 
individual  reactions  to  be  in  equilibrium  or  some  of  the  species  concentration  can  be 
assumed  to  be  in  quasi-steady  state  [5]. 


G.  EQUIVALENCE  RATIO  OF  FUEL-OXIDIZER  MIXTURE 

Typically,  combustion  of  fuel-oxidizer  mixtures  is  specified  in  terms  of  the 
equivalence  ratio,  cj).  The  equivalence  ratio  is  independent  of  the  mass,  m  or  the  number 
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of  moles,  n  that  is  used  to  eharaeterize  the  proportion  of  the  fuel  to  oxidizer  mixture.  It  is 
defined  as  the  ratio  of  the  fuel-oxidizer  ratio  to  the  stoiehiometrie  fuel-oxidizer  ratio  and 
ean  be  mathematieally  represented  as  follows: 


The  eomposition  and  thermodynamie  properties  of  fuel-oxidizer  mixtures  during 
eombustion/detonation  are  dependent  on  the  equivalenee  ratio  and  eonsequently  the 
detonation  eell  width  if  the  wave  propagates  at  supersonie  speeds.  A  ratio  greater  than  1 
would  imply  exeess  fuel  or  fuel  rieh  mixtures  than  would  be  required  for  a  stoiehiometrie 
reaetion  and  eonversely  a  ratio  less  than  1  would  imply  exeess  oxidizer  or  fuel  lean 
mixtures. 

H.  COMPUTATIONAL  NUMERICAL  SCHEMES 

In  eomputational  fluid  dynamies,  there  are  primarily  two  main  sehemes,  explieit 
and  implieit  sehemes  implemented  in  order  to  diseretize  the  resulting  differential 
equations  and  resolve  the  required  flow  field.  Explieit  or  forward  differeneing  sehemes 
suffer  from  eonditional  stability  issues  where  a  solution  may  blow  up  if  the  speeified  time 
step  is  too  large  and/or  the  eomputational  eell  size  is  too  small.  Implieit  sehemes  are 
inherently  stable  as  the  forward  solution  is  eoupled  together  with  the  previous  solution 
and  henee  will  not  generate  an  output  that  will  inerease  to  infinitely  large  values. 
Sehemes  of  sueh  types  will  be  employed  as  the  main  seheme  in  study  for  simulation  of 
detonation  waves  in  PDE  eonfigurations.  Eigure  1 1  show  a  eomputational  moleeule  based 
on  an  implieit  seheme  where  p  defines  the  spatial  grid  points  in  spaee  and  q  defines  the 
temporal  grid  points  on  an  x-t  diagram. 
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For  a  mathematical  representation  of  an  implicit  scheme,  consider  a  simple  one¬ 
dimensional  scalar  equation: 


du  du 

—  +  a —  =  0 
dx  dt 


(14) 


Discretization  of  the  above  equation  based  on  an  implicit  scheme  results  in  the  following: 


(15) 


- is  also  known  as  the  Courant-Fredrichs-Lewy  or  CFL  number.  Numerical  solution  of 

Ax 

the  resulting  flow  field  can  be  specified  either  in  terms  of  a  fixed  time  step  or  the  CFL 
number  where  the  time  step  is  computed  based  on  the  computational  mesh  size.  For  a 
sufficiently  small  fixed  time  step,  specification  of  the  CFL  number  will  result  in  faster 
numerical  computations  [6]. 
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III.  SIMULATION  SETUP  (GRID  SENSITIVITY) 


A,  WAVE  IMPLOSION 

Wave  implosion  is  a  technique  that  focuses  a  shock  front  to  generate  a  region  of 
high  pressure  and  temperature  gases  due  to  adiabatic  compression  of  gas  products  as  it 
flows  into  a  region  of  decreasing  area  [1],  As  a  result,  extremely  high  energy  density 
localized  regions  are  formed  due  to  interacting  compression  waves  and  this  raises  the 
pressure  of  a  typical  fuel-oxidizer  mixture  above  the  C-J  pressure,  resulting  in  an 
overdriven  detonation  wave.  Figure  12  illustrates  an  example  of  a  collapsing  toroidal 
wave  front  implosion.  This  concept  is  being  explored  to  increase  specific  thrust  of  PDFs 
by  harnessing  the  extremely  high  pressures  as  a  result  of  the  overdriven  imploding  shock 
front  to  increase  the  maximum  diameter  for  a  DDT  process  assuming  a  reasonable  run 
distance  required  for  transition  to  detonation.  A  commercial  Computational  Fluid 
Dynamics  (CFD)  platform  used  to  perform  simulation  studies  on  the  implosion  concepts 
for  PDFs  is  CFD-l-l-  from  Metacomp  Technologies.  It  is  widely  used  by  commercial 
companies  such  as  General  Flectric  and  Air  Force  Research  Faboratories  (AFRF)  as  a 
computational  tool  for  development  on  proof  of  concept  PDFs  and  also  supported  by 
Office  of  Naval  Research  (ONR). 


PI 

P2 

P3 

P4 


Figure  12,  Collapsible  toroidal  wave  implosion  [From  1] 
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B.  HIGH  FIDELITY  SIMULATIONS  (GRID  SENSITIVITY  ANALYSIS) 


In  order  to  explore  implosion  concepts,  high  fidelity  computational  fluid  dynamic 
(CFD)  models  based  on  reduced  chemical  kinetic  mechanisms  are  pertinent  for  accurate 
and  robust  simulation  to  describe  and  analyze  the  physics  of  detonation  wave  propagation 
in  pulse  detonation  engines.  High  fidelity  models  are  sensitive  to  the  induction  length  for 
a  specific  fuel-oxidizer  mixture  such  that  the  transient  flow  kinetics  and  shock  physics 
can  be  accurately  captured  as  shown  by  the  ZND  model  in  the  previous  section. 

Using  a  test  case  model  for  Hydrogen-Air  (H2-Air)  mixtures,  a  grid  sensitivity 
analysis  was  performed  on  a  typical  PDE  tube  configuration  using  a  kinetic  reaction 
model  based  on  a  reduced  set  of  18  reactions  and  9  species  [7]  and  [8]  highlighted  in 
Table  2. 


Table  2,  Reduced  chemical  kinetic  mechanism  for  H2-Air  mixture 


Number 

Reaction 

A" 

n" 

F  ^ 

1 

H2  +  O2  ^  20H 

0.17ell 

0 

0.2015e9 

2 

O2  +  H  ^  20H 

0.142el2 

0 

0.6866e8 

3 

H2  +  OH  ^  H2O  +  H 

0.316e5 

1.8 

0.12686e8 

4 

H2  +  0  ^  OH  +  H 

0.207el2 

0 

0.57568e8 

5 

20H  ^  H2O  +  0 

0.55ell 

0 

0.29307e8 

6 

OH  +  H  +  M  ^  H2O  +  M 

0.2211el7 

-2 

0 

7 

H  +  H  +  M  ^  H2 

0.653el2 

-1 

0 

8 

O2  +  H  +  M  ^  HO2  +  M 

0.32el3 

-1 

0 

9 

OH  +  HO2  ^02  +  H2O 

0.5ell 

0 

0.4186e7 

10 

H  +  HO2  ^  H2  +  O2 

0.253ell 

0 

0.29307e7 

11 

H  +  HO2  ^  20H 

0.199el2 

0 

0.7536e7 

12 

0  +  HO2  ^  O2  +  OH 

0.5ell 

0 

0.4186e7 

13 

2HO2  ^  O2  +  H2O2 

0.199el0 

0 

0 

14 

H2  +  HO2  ^  H  +  H2O2 

0.301e9 

0 

0.7829e8 

15 

OH  +  H2O2  ^  H  +  H2O2 

0.102ell 

0 

0.7954e7 

16 

H  +  H2O2  ^  H2O  +  OH 

0.5el2 

0 

0.4186e8 

17 

0  +  H2O2  ^  OH  +  HO2 

0.199ell 

0 

0.247e8 

18 

H2O2  +  M  ^  20H  +  M 

0.121el5 

0 

0.19049e9 

in  units  of  kmofm^,  s,  K  and  JAmol 
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Nitrogen  and  carbon  chemistry  reactions  were  excluded  in  the  reaction  model  so 
as  to  optimize  the  required  computational  effort.  Only  slight  differences  in  CJ  detonation 
wave  speeds  and  temperatures  for  typical  hydrocarbon  fuel-air  mixtures  on  the  order  of 
6%  to  9%  respectively  [9]  were  observed  as  compared  to  the  full  kinetic  model 
comprising  the  full  nitrogen  and  carbon  chemistry  for  Jet  A/ Air  mixture,  as  shown  in 
Figures  13  and  14,  which  was  deemed  as  acceptable. 
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Figure  13,  Comparison  of  detonation  velocities  for  Jet  A/ Air  mixture  [From  9] 
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Comparison  of  C-J  temperatures  for  Jet  A/ Air  mixtures  [From  9] 
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Figure  14, 


c. 


DETONATION  CELL  WIDTH  OF  H2-AIR  FOR  DIFFERENT 
EQUIVALENCE  RATIOS 


Determination  of  the  eomputational  mesh  size  for  the  model  was  foeused  mainly 
on  stoiehiometrie  Hi-Air  mixtures.  This  is  beeause  it  is  well  defined  in  literature  that 
experimental  results  at  whieh  detonation  eell  width  oeeurs  at  a  minimum  eorresponds  to 
=  1  as  indieated  in  Figure  15  refereneed  from  [10],  [11],  [12]  and  [13]  for  a  H2-Air 
mixture.  The  induetion  length  is  in  turn  proportional  to  the  detonation  cell  width  and  the 
relation  is  typically  a  linear  one  with  a  proportionality  constant  for  different  fuel-oxidizer 
mixtures.  Therefore,  this  would  present  the  baseline  mesh  size  required  for  accurate 
simulation  of  detonation  shock  physics. 


Qril  aWlk  ucVk  IR-Ak 


Figure  15,  Graph  of  detonation  cell  width  vs  equivalence  ratio  for  H2-Air 

mixture  [From  14] 
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D, 


PDE  CONFIGURATION 


For  the  grid  sensitivity  analysis,  a  2-D  axisymmetrie  eylinder  model  was 
generated  on  the  Solid  Works  platform  with  the  dimensions  outlined  in  Figure  16.  A 
detonation  wave  eonsisting  of  the  detonation  produets  was  assumed  to  be  initiated  at  CJ 
Pressure  and  CJ  Temperature  eonditions  eharaeteristie  of  the  fuel-oxidizer  mixture  as  part 
of  the  initial  eonditions.  Mass  fraetion  of  the  detonation  produets  are  obtained  from 
NASA  Chemieal  Equilibrium  Analysis  (CEA)  [15]  and  ineluded  in  Appendix  A.  The 
reaetant  initial  eonditions  were  assumed  to  be  at  1  atmosphere  and  298  K,  whieh  would 
mimie  experimental  eonditions  for  further  validation. 


Figure  16,  2-dimensional  axisymmetrie  cylinder  model 


The  model  was  gridded  on  the  Multipurpose  Intelligent  Meshing  Environment 
(MIME)  platform  for  an  unstruetured  mesh  (triangles)  topology,  whieh  was  optimized  as 
a  pre-proeessor  for  CED-H-I-.  It  would  have  been  ideal  to  resolve  the  flow  field  in  terms  of 
a  struetured  mesh  topology  (quadrilaterals)  so  as  to  minimize  the  eomputational 
workload;  however  MIME  did  not  have  the  eapability  to  generate  struetured  grids.  An 
illustration  of  the  cylinder  model  for  unstructured  mesh  topology  can  be  found  in  Eigure 
17. 


Length  (m) 
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Figure  17.  2-D  axisymmetric  model  (Unstructured  mesh  size) 


For  the  implosion  study,  a  modified  PDF  eonfiguration  was  eoneeived  based  on  a 
2-D  eylinder  model  as  depleted  in  Figure  18  with  an  axi-symmetrie  axis.  This  model 
utilizes  the  eoneept  of  implosion  with  an  assumption  of  a  starting  detonation  wave 
propagating  at  the  inlet.  The  lead  shoek  front  of  the  detonation  wave  should  be  sustained 
by  the  implosion  eoneept  and  will  not  deeay  to  a  rarefaetion  wave  as  it  expands  into  the 
tube  spaee.  Further  parametrie  studies  are  then  eondueted  by  inereasing  the  radius  of  the 
model  until  the  detonation  wave  fails,  thereby  eompromising  the  DDT  length  that  is 
required  for  sueh  geometrie  eonfigurations. 


blat  CPnduct^ 


E,  SIMULATION  PARAMETERS 

The  simulation  parameters  used  to  perform  the  grid  sensitivity  analysis  is 
summarized  in  Table  3.  The  unsteady  eompressible  Euler  equation  (no  viseosity  terms) 
was  ehosen  as  the  main  equation  set  governing  the  numerieal  eomputations.  The 
parameter  of  interest  was  to  evaluate  the  thermodynamie  parameters  and  wave  speed 
along  the  symmetry  axis  so  as  to  eapture  the  detonation  wave  physies  rather  than  near  the 
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walls  (where  viscous  effects  become  significant).  A  Total  Variation  Diminishing  (TVD) 
limiter  scheme  based  on  the  minmod  function  was  used  in  order  to  prevent  spurious 
oscillation  of  the  resolved  flow  field  in  the  vicinity  of  strong  gradients,  which  represent 
shocks.  Near  flow  field  regions  where  shocks  are  encountered,  they  switch  the  spatial 
discretization  scheme  to  a  first  order  accurate  method  while  second-order  accurate 
numerical  schemes  are  computed  for  the  rest  of  the  flow  field. 


Table  3,  Simulation  setup  on  CFD-i-i-  platform 


Equation  Set  Type 

i) 

Unsteady  Compressible  Euler 

ii) 

Equation  of  State:  Ideal  Gas 

Riemann  Solver 

i) 

Minimum  Dissipation:  EHS  &  RHS 

ii) 

Activate  Pressure  Switch:  Supersonic 

iii) 

Activate  Pressure  Gradient  Normal 

Boundary  Conditions  (See  Figure  17) 

i) 

Top:  Multi-species  inviscid  surface  tangency  (Wall) 

ii) 

Eeft:  Multi-species  inviscid  surface  tangency  (Wall) 

iii) 

Right:  Multi-species  inviscid  surface  tangency  (Wall) 

iv) 

Bottom:  X-axis  symmetry 

Boundary  Conditions  (See  Figure  18) 

i) 

Inlet:  Specified  by  inlet  stagnation  temperature,  pressure  and 
detonation  products 

ii) 

Wall:  Multi-species  inviscid  surface  tangency 

iii) 

Symmetry:  X-axis  symmetry 

Time  Integration 

i)  Implicit 

Spatial  Discretization 

i) 

2““*  Order  Accuracy  in  Space 

ii) 

Dimension  of  polynomial:  2-D  axisymmetric 

iii) 

Axis  stability  enhancement:  Yes 

iv) 

Types  of  Total  Variation  Diminishing  (TVD)  Eimiter:  Minmod 
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IV.  DISCUSSION  OF  RESULTS  AND  ANALYSIS 


A.  INITIAL  COMPARISON  OF  DIFFERENT  GRID  SIZES 

A  grid  size  from  0.25  mm  to  1.0  mm  was  first  considered  in  the  grid  sensitivity 
analysis  so  as  to  determine  a  computational  mesh  size  that  is  sufficiently  small  enough  to 
capture  the  physics  of  the  detonation  wave.  The  ZND  induction  length  for  detonation  of 
stoichometric  H2-Air  mixtures  (0.23  mm)  was  used  as  a  baseline  comparison  [16]. 

As  outlined  previously  in  the  PDE  configuration  description,  a  detonation  wave 
was  assumed  to  propagate  at  the  start  with  conditions  pertaining  to  C-J  pressure  and 
temperature  obtained  from  NASA  CEA  code.  This  is  outlined  in  Table  4.  (Note  that  the 
initial  reactants  were  assumed  to  be  at  Pressure  =  1  atm  and  Temperature  =  298K).  This 
would  mimic  experimental  conditions  for  purpose  of  validation.  The  single  step  reaction 
equation  is  indicated  in  Table  5. 

Table  4,  C-J  parameters  for  stoichometric  H2-Air  mixtures 


Vcj 

(m/s) 

Tcj  (K) 

Pcj  (MPa) 

1.0 

1964.8 

2942 

1.58 

Table  5,  Single-step  reaction  for  stoichometric  Hi-Air  mixtures 


c|j 

Reaction 

1.0 

H2  +  14(02  +  3.76N2)  ->  H2O  +  3.76/2N2 

An  implicit  scheme  with  a  fixed  time  step  of  1  nanosecond  was  first  used  during 
the  initial  comparison  process.  In  order  to  validate  the  accuracy  of  the  simulation  run  for 
each  computational  mesh  size,  the  computed  detonation  velocity  between  each  time  step 
along  the  symmetry  axis  was  computed  using  a  first  order  finite  difference  method  and 
compared  against  the  C-J  velocity.  The  C-J  velocity  is  the  velocity  at  which  a  steady  state 
detonation  wave  propagates,  which  has  been  verified  in  shock  tube  experiments,  as  noted 
by  Pickett  and  Davis  [17]. 
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1,  Computation  of  Detonation  Velocity 

Detonation  velocity  was  determined  by  noting  the  location  of  the  shock  front  in 
the  form  of  a  pressure  discontinuity  along  the  symmetry  axis  of  the  cylinder 
configuration  between  each  subsequent  time  step  as  illustrated  in  Figures  19  and  20. 


Figure  19,  Pressure  contour  for  H2-Air  mixture  at  3  microsecond  time  step 


Figure  20,  Pressure  contour  for  Hi-Air  mixture  at  4  microsecond  time  step 
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To  obtain  the  wave  speed,  a  five-point  steneil  first  order  finite  differenee  method 
[18]  as  shown  in  Figure  21  is  performed  in  one-dimension  where  x  represents  the  spatial 
diseretized  loeation  of  the  shoek  front  and  t  represents  the  fixed  time  step  pertaining  to 
the  shoek  front  loeation. 


O 


“O 

K(H-2) 


Figure  21,  Schematic  of  the  5-point  stencil  method 


The  detonation  wave  speed,  is  then  obtained  via: 


V, 


DET 


-  x{t  +  7)+  8x(t  + 1)  -  8x(t  - 1)  +  x{t  -  2) 

uKt 


(16) 


2,  Results  from  Initial  Comparison 

Results  of  the  normalized  detonation  veloeity  with  C-J  veloeity  based  on  eell 
mesh  sizes  of  0.25  mm  (~  One  ZND  length),  0.5  mm  (~  Two  ZND  lengths)  and  1.0  mm 
(~  Four  ZND  lengths)  for  stoiehiometrie  H2-Air  mixtures  against  run  distanee  along  the 
symmetry  axis  of  the  PDF  cylinder  is  plotted  in  Figure  22.  Convergence  plots  for  each  of 
the  simulation  runs  pertaining  to  each  individual  mesh  size  are  included  in  Appendix  B  to 
verify  that  the  resolved  flow  field  solutions  are  bounded. 
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Plot  of  Wave  Velocitv/CJ  Velocity  against  Run  Distance  for  Different 

Mesh  Size 
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Figure  22.  Plot  of  wave  velocity/CJ  velocity  against  run  distance  for  different 

mesh  sizes 


All  3  cell  sizes  indieate  good  agreement  with  the  predicted  C-J  veloeity  towards 
the  end  of  the  PDE  cylinder.  The  detonation  waves  decay  towards  the  steady  state  C-J 
velocity  at  the  end  where  the  order  of  aecuracy  from  the  mean  C-J  velocity  ranged  from 
3%  to  5%.  Table  6  gives  a  statistical  summary  of  the  detonation  wave  velocities  in  terms 
of  the  standard  deviation  and  the  average  normalized  detonation  wave  speeds.  This 
elearly  implies  that  for  simulation  of  a  fully  developed  detonation  wave  propagating  at  C- 
J  velocity,  the  cell  sizes  are  clearly  sufficient.  However,  the  speeified  initial  conditions 
were  not  enough  to  induce  the  ‘wave’  transient  two-dimensional  effects  constituting  the 
longitudinal  and  transverse  wave  interactions  that  have  been  observed  in  PDE  detonation 
soot  foil  imprint  experiments.  Hence,  further  grid  sensitivity  studies  were  performed  with 
the  initial  detonation  wave  propagating  at  temperature  and  pressure  conditions  specified 
by  Zaev  et  al  [8]. 
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Table  6, 


Set  of  evaluation  parameters  for  different  mesh  sizes 


Mesh  Size 
(mm) 

Average 
Detonation 
Wave  Velocity 
(m/s) 

Standard  Deviation 
of  Detonation  Wave 
Velocity  (m/s) 

Average 

Detonation  Wave 
Velocity  /  C  J 
Velocity 

0.25  (1  ZND  Length) 

1892.58 

70.88 

0.96 

0.50  (2  ZND  Lengths) 

1864.54 

77.87 

0.95 

1.00  (4  ZND  Lengths) 

1898.81 

175.55 

0.97 

Also,  the  computational  effort  required  to  complete  a  simulation  run  for  0.25  mm 
took  8  days  based  on  an  implicit  scheme  with  a  fixed  time  step.  To  reduce  and  optimize 
the  computational  time,  an  implicit  scheme  with  CFL  number  equal  to  one  was  adopted. 
In  order  to  validate  the  implicit  scheme  based  on  CFL  number,  the  same  test  case  model 
using  a  mesh  size  of  0.25  mm  was  simulated  using  these  two  employed  schemes  and 
shown  in  Figure  23.  The  results  indicate  good  agreement  with  each  other  with  the 
numerical  accuracy  not  compromised  by  adopting  a  different  scheme. 


Figure  23,  Plot  of  wave  velocity/CJ  velocity  against  run  distance  for  different 

implicit  time  schemes 
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B. 


FURTHER  GRID  SENSITIVITY  ANALYSIS 


Further  analysis  was  performed  on  cell  sizes  ranging  from  0.25  mm  (~  One  ZND 
length),  0.0625  mm  (~  One-quarter  ZND  length)  and  0.05  mm  (~  One-fifth  ZND  length). 
A  plot  of  the  normalized  wave  speeds  along  the  symmetry  axis  with  the  C-J  velocity  as  a 
function  of  the  run  distance  of  the  cylinder  is  shown  in  Figure  24.  Initial  conditions  for  a 
starting  detonation  wave  are  specified  at  a  pressure  of  3.95  MPa  and  a  temperature  of 
2942  K. 


1.60 


Plot  of  Wave  Velocitv/CJ  Velocity  against  Run  Distance  for  Different 

Mesh  Size 


- Wave  Velocity/CJ  Velocity  (Mesh  Size  0.125  mm) 


Run  Distance(m) 


Figure  24,  Plot  of  wave  velocity/CJ  velocity  against  run  distance  for  different 

mesh  sizes  for  Courant  number  of  1 


At  a  computational  mesh  size  of  0.0625  mm,  initial  wave  transients  could  be 
observed  for  an  under-driven  wave  propagating  at  a  sub-CJ  speed  and  increase  after  a 
small  run  distance  to  an  over-driven  detonation  wave  before  decaying  and  oscillating 
about  the  steady  state  C-J  velocity.  This  clearly  shows  the  transition  to  detonation  or  the 
DDT  process.  Decreasing  the  cell  size  from  0.0625  mm  to  0.05  mm  shows  no  change  in 

the  initial  transients  as  well  as  the  decay  behavior  to  the  steady  state  C-J  velocity,  which 

32 


allows  us  to  determine  that  0.0625  mm  would  be  a  praetieal  eell  size  to  eapture  the 
detonation  physios,  as  well  as  apply  to  other  fuel-oxidizer  mixtures. 

Further  insight  oould  also  be  gained  by  first  analyzing  the  pressure  and 
temperature  oontours  as  well  as  the  histories  for  the  given  mesh  size  of  0.0625  mm  at  the 
mid-point  of  the  oylinder  oonliguration  as  indioated  in  Figures  25,  26,  27  and  28 
respeotively. 

Behind  the  shook  front,  the  detonation  wave  struoture  oomprising  of  transverse 
wave  interaotions  along  with  the  longitudinal  propagating  wave  front  are  observed,  whioh 
are  olearly  multi-dimensional  effeots  for  both  pressure  and  temperature  oontours.  The  thin 
shook  front  layer  is  propagating  at  a  pressure  state  of  2.4  MPa  that  is  higher  than  the  C-J 
pressure  of  1.58  MPa.  Just  behind  the  shook  front,  the  detonation  wave  struoture 
approximately  45  mm  long  oomprising  of  the  triple  points  are  olose  ooupled,  whioh  is 
fueling  the  detonation  wave  prooess  followed  by  an  expansion  zone.  Similarly,  the  H2-Air 
mixture  reaots  at  an  almost  oonstant  C-J  temperature  oondition  of  3000  K  throughout  the 
detonation  wave  struoture  zone. 

From  the  pressure  history  data,  the  pressure  state  for  lead  shook  front  relaxes  to 
1.0  MPa  before  inoreasing  baok  up  to  2.0  MPa.  The  initial  deorease  is  due  to  the 
interaotion  of  the  transverse  wave  propagating  from  the  side  walls  and  the  symmetry  axis 
of  the  oylinder.  However,  when  the  triple  points  interseot,  the  waves  reinforoe  eaoh  other 
and  the  superposition  prooess  inoreases  the  looalized  pressure  points.  This  promotes 
oonditions  that  sustain  the  propagation  of  the  detonation  wave. 
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Figure  25.  H2  -  Air  pressure  contour  plot  at  65  mm 


Figure  26.  H2-  Air  pressure  history  plot  at  65  mm 
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Figure  27.  H2  -  Air  temperature  contour  plot  at  65  mm 
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Figure  28.  H2  -  Air  temperature  history  plot  at  65  mm 
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The  detonation  wave  is  sustained  throughout  by  analyzing  the  contour  and 
histories  pertaining  to  conditions  of  pressure  and  temperature  at  the  cylinder  end.  This  is 
illustrated  in  Figures  29,  30,  31  and  32.  The  same  detonation  wave  structure  closely 
coupled  to  the  shock  front  is  observed  with  no  change  in  the  reaction  zone  length  with  a 
larger  expansion  zone.  The  lead  shock  front  is  propagating  at  a  slightly  lower  pressure 
state  of  2.3  MPa.  As  before,  the  pressure  decreases  before  increasing  back  up  to  2.4  MPa 
due  to  the  interaction  of  the  wave  transients.  The  reaction  zone  conditions  are 
approximately  constant  at  3000  K. 

One  interesting  feature  to  take  note  is  when  the  detonation  wave  structure  is 
closely  analyzed,  as  shown  in  Figure  33,  the  detonation  cell  width,  which  characterizes 
the  triple  point  intersection  of  the  mach  stem,  incident  shock  and  reflected  shock  is 
approximately  8  mm,  which  was  taken  to  be  the  largest  cell  size  since  the  detonation  cell 
width  is  irregular  in  nature.  This  coincides  with  experimental  observations  as  predicted  in 
Figure  15. 
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Figure  29,  H2  -  Air  pressure  contour  plot  at  122  mm 
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Figure  30,  Hi-Air  pressure  history  plot  at  122  mm 
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Figure  31.  H2  -  Air  temperature  contour  plot  at  122  mm 


37 


Symmetry  Length  (m) 


Figure  32,  H2  -  Air  temperature  history  plot  at  122  mm 


Figure  33,  Close-up  view  of  detonation  wave  structure 
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This  concludes  the  grid  sensitivity  analysis  and  establishes  that  the  eell  size  of 
0.0625  mm  (-'  one-quarter  induetion  length)  is  suffieient  to  eapture  the  detonation  wave 
physies  in  a  PDE  eylinder  eonfiguration.  The  same  eomputational  eell  size  will  be  used  in 
subsequent  simulations  with  a  modified  eonfiguration  for  analysis  of  axi-symmetrie 
toroidal  implosions. 

C.  IMPLOSION  ANALYSIS 

Using  the  modified  PDE  eonfiguration  illustrated  in  Eigure  18,  a  detonation  wave 
was  assumed  to  propagate  at  the  inlet  with  stagnation  pressure  and  temperature  equal  to 
3.95  MPa  and  2942  K  with  the  detonation  produets  obtained  from  NASA  CEA.  The 
reaetant  eonditions  were  speeified  to  be  at  an  initial  pressure  of  1  atm  and  temperature  of 
298  K.  This  was  similar  to  the  eonditions  simulated  for  the  grid  sensitivity  analysis.  The 
eell  size  used  was  0.0625  mm  as  determined  previously,  whieh  was  required  to  eapture 
the  eomplex  flow  field  governing  the  shoek  physies.  The  same  test  ease  model  for 
stoiehiometrie  reaetion  of  H2-Air  mixture  assuming  a  ehemieal  kinetie  reaetion  setup  of 
18  reaetions  and  9  speeies  listed  in  Table  2  was  used  for  the  implosion  analysis. 

In  order  to  determine  whether  the  propagating  detonation  wave  ean  be  sustained 
by  the  implosion  proeess,  the  detonation  veloeity  was  evaluated  along  the  mid-axis  of  the 
inlet  (whieh  would  be  the  eenter  of  the  implosion  event)  and  along  the  symmetry  axis  of 
the  shoek  tube  eonfiguration  as  shown  in  Eigure  34.  The  resultant  detonation  veloeity  was 
eompared  against  the  C-J  veloeity  to  see  if  the  wave  ean  persist  as  a  detonation  wave. 
Therefore,  the  eorresponding  analysis  are  deeoupled  into  two  portions,  the  first  pertaining 
to  the  mid-axis  inlet  (implosion)  and  the  seeond  pertaining  to  the  symmetry  axis. 
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Figure  34,  Evaluation  of  C-J  wave  speed  for  implosion  study 


1,  Implosion  Along  the  Mid-axis  Inlet 

Results  of  the  normalized  detonation  wave  veloeity  to  C-J  veloeity  versus  run 

distanee  along  the  mid-axis  of  the  inlet  were  plotted  in  Figure  35.  From  the  graph,  the 

detonation  wave  veloeity  deereases  initially  up  to  a  run  distanee  of  15  mm  (-'  0.4  of  the 

total  radius)  below  C-J  veloeity  as  it  travels  along  the  radial  direetion  towards  the 

symmetry  axis.  This  ean  be  attributed  to  rarefaetion  effeets  as  the  detonation  wave  starts 

to  expand  from  the  inlet  entranee.  This  ean  explieitly  illustrated  from  the  respeetive 

snapshot  pressure  and  temperature  eontours  in  Figures  36  and  37  for  the  imploding  wave 

at  a  run  distanee  of  10  mm.  From  Figure  36,  there  is  an  expansion  zone  whereby  the 

starting  detonation  wave  deereases  in  pressure  from  speeified  inlet  stagnation  eonditions 

along  the  inlet  axis  all  the  way  to  the  reaetion  zone  attaehed  to  the  lead  shoek  front  that  is 

propagating  at  the  wave  veloeity.  This  reaetion  zone  is  a  region  of  eompression  where  the 

pressures  start  to  build  up  fueled  by  the  ehemieal  reaetion  of  the  fuel-oxidizer  mixture. 

Notiee  that  this  zone  is  loealized  only  at  the  tip  of  the  spherieal  wave  front  due  to  the 

nature  of  the  toroidal  imploding  geometry  and  is  a  typieal  example  of  a  shook  foousing 

teohnique.  There  is  a  large  high  temperature  region  approximately  equal  to  C-J 

temperature  throughout  the  spherieal  shook  wave  exoept  for  the  region  immediately  to  the 

left  and  right  of  the  inlet  entranee  whereby  no  oommunioation  exists  between  the  flow 

field  within  the  shook  and  reaotant  eonditions.  The  observed  high  temperature  region  is 

supported  by  the  heat  release  due  to  the  reaetion  of  the  fuel-oxidizer  mixture  due  to  the 
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interaction  of  the  multi-dimensional  wave  structure  as  the  shock  front  implodes  towards 
the  symmetry  axis.  The  pressures,  temperatures  as  well  as  the  kinetic  chemical  reaction 
rates  are  clearly  driving  the  spherical  wave  front. 


Plot  of  Wave  Velocitv/CJ  Velocity  against  Run  Distance  for  Inlet  Axis 


Tube  Radius  (m) 


Figure  35,  Plot  of  wave  velocity/CJ  velocity  against  tube  radius  along  mid-axis  of 

inlet 


41 


0.1 


CO 


CO 

CC 


P 

3,00E  +  06 
2,71  E  +  06 
2,42E  +  06 

2,1  3E  +  06 

1  .84E  +  06 

1  .55E  +  06 

1  .26E  +  06 

0.08 

1 

0.06 

LJ  ,  O  LJ  +  LI  o 

3,90E  +  05 

1  ,00E-i-05 

0.04 


0.02 


°0  0.02  0.04  0.06  0.08  0.1 

Symmetry  Length  (m) 


Figure  36,  Pressure  contour  for  imploding  geometry  at  run  distance  of  10  mm 

along  inlet  axis 
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Figure  37,  Temperature  contour  for  imploding  geometry  at  run  distance  of  10 

mm  along  inlet  axis 
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As  the  wave  propagates  beyond  15  mm  towards  the  symmetry  axis,  the  wave 
velocity  starts  to  increase  back  upwards  towards  C-J  velocity  by  virtue  of  the  imploding 
geometry.  This  illustrates  the  fact  that  the  detonation  wave  is  continuously  supported  by 
the  shock  focusing  process  that  can  be  determined  from  the  C-J  velocity.  The  pressure 
and  temperature  contours  corresponding  to  a  run  distance  of  20  mm  are  illustrated  in 
Figures  38  and  39  respectively.  The  spherical  wave  front  from  the  previous  run  distance 
of  10  mm  has  now  transformed  into  little  spherical  wavelets  coalescing  into  the  leading 
edge  of  the  detonation  wave.  This  is  due  to  little  compression  waves  that  are  accelerated 
due  to  generation  of  heat  release  from  the  reaction  of  the  fuel-oxidizer  mixture.  This  heat 
release  form  high  temperature  regions  that  increases  the  sonic  velocity  of  the  gas 
products  tries  to  catch  up  with  the  leading  edge  of  the  shock  front.  Also,  extremely  high 
localized  pressure  regions  denoted  by  tiny  red  regions  are  also  observed  at  the  leading 
edge  that  is  at  least  a  factor  of  two  greater  than  the  C-J  pressure.  A  look  at  the 
temperature  contour  reveals  an  interesting  phenomenon,  where  transverse  wave 
interactions  are  observed  from  these  spherical  wavelets  at  the  protruding  edge  of  the  lead 
shock  front  extending  for  a  distance  of  10  mm  along  the  mid-axis  inlet  as  shown  in  the 
close-up  view  for  the  same  temperature  contour  in  Figure  40.  There  is  also  an  extremely 
thin  region  separating  the  ambient  temperature  from  the  high  temperature  region  - 
characteristic  of  a  shock  leading  front.  Further  away  from  the  mid-axis  of  the  inlet  on 
both  sides,  the  thickness  of  this  region  increases,  which  indicates  that  the  reactions  are 
falling  behind  the  shock  front,  and  the  detonation  process  is  decaying  due  to  expansion 
effects. 
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Figure  38,  Pressure  contour  for  imploding  geometry  at  run  distance  of  20  mm 

along  inlet  axis 
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Figure  39,  Temperature  contour  for  imploding  geometry  at  run  distance  of  20 

mm  along  inlet  axis 
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Figure  40,  Close-up  view  of  temperature  contour 

The  pressure  and  temperature  history  plots  along  the  mid-axis  inlet  are  used  to 
quantify  the  high  pressures  and  temperatures  driving  the  imploding  wave  as  shown  in 
Figures  41  and  42.  The  lead  shock  front  is  propagating  at  a  pressure  of  2.5  MPa,  which 
essentially  represents  an  overdriven  detonation  wave.  Immediately  behind  the  shock 
front,  the  pressure  drops  drastically  to  1.2  MPa  followed  by  a  constant  pressure  region 
approximately  8  mm  thick.  As  for  the  temperature  at  the  leading  edge,  it  increases  slowly 
from  2500  K  in  an  oscillating  manner  to  almost  2900  K,  which  equals  approximately  to 
CJ  temperature  before  spiraling  downwards  to  2000  K.  This  spiraling  behavior  is  due  to  a 
huge  expansion  effect  where  the  pressure  decreases  from  the  constant  pressure  region  to 
ambient  pressure  conditions. 
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Figure  41,  Pressure  history  of  imploding  wave  at  run  distance  of  20  mm  along 

inlet  axis 
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Figure  42,  Temperature  history  of  imploding  wave  at  run  distance  of  20  mm 

along  inlet  axis 
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As  the  wave  eontinues  to  implode  towards  the  symmetry  axis,  the  loealized  high 
pressure  region  attaehed  to  the  shoek  front  expand  radially  outwards  with  the  wave  at  a 
run  distanee  of  38  mm,  which  is  almost  the  radius  of  the  PDE  tube  configuration  as 
shown  in  Figure  43.  The  pressure  history  in  Figure  44  verifies  the  fact  that  the  pressure 
continuously  increases  beyond  the  C-J  pressure  to  3.5  MPa  and  this  further  reinforces  and 
builds  up  the  detonation  wave  speed.  Similarly,  the  expansion  wave  region  thickness 
increases  as  the  high  pressure  relieves  itself  to  0.9  MPa  before  a  sharp  decrease  to 
ambient  pressure  conditions.  This  indicates  that  a  second  shock  front  is  formed  in  front 
of  a  set  of  supersonic  jets  of  products  leaving  the  inlet  region  and  entering  the  detonation 
tube.  This  phenomenon  will  be  treated  robustly  in  future  endeavors. 
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Figure  43,  Pressure  contour  for  imploding  geometry  at  run  distance  of  38  mm 

along  inlet  axis 
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Figure  44,  Pressure  history  of  imploding  wave  at  run  distance  of  38  mm  along 

inlet  axis 

Similar  to  the  observed  pressures,  there  is  a  radial  expansion  of  the  high 
temperature  region  eorresponding  to  approximately  C-J  temperature  outwards  from  the 
inlet  mid-axis  as  illustrated  in  the  pressure  eontour  of  Figure  45.  The  temperature  history 
plot  shows  an  almost  eonstant  temperature  region  of  2900  K  extending  from  the  shoek 
front  to  a  thiekness  of  20  mm,  whieh  is  almost  half  the  radius  of  the  PDF  tube. 

Analyzing  the  entire  propagation  of  the  toroidal  wave  front  motion  proves  that 
even  though  the  detonation  wave  initially  loses  its  speed  due  to  sudden  expansion  effeets, 
the  implosion  proeess  eoupled  with  the  reaetion  of  the  fuel-oxidizer  mixture  strengthens 
and  inereases  the  wave  veloeity  as  it  propagates  towards  the  symmetry  axis.  The  pressure 
state  just  before  the  symmetry  also  shows  that  the  initial  wave  builds  up  to  an  overdriven 
detonation  shook  wave  propagating  by  a  faotor  of  2.2  of  the  C-J  pressure  due  to  the 
implosion  proeess. 
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Figure  45,  Temperature  contour  for  imploding  geometry  at  run  distance  of  38 

mm  along  inlet  axis 


Figure  46,  Temperature  history  for  imploding  geometry  at  run  distance  of  38 

mm  along  inlet  axis 
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2,  Implosion  Analysis  Along  Symmetry  Axis 

Looking  at  the  passage  of  the  shock  front  along  the  symmetry  axis,  it  starts  out  at 
a  highly  overdriven  state  with  a  wave  velocity  propagating  by  a  factor  of  1.8  times  the  C- 
J  velocity.  This  is  a  highly  transient  state  that  will  eventually  decay  to  the  steady  state  C-J 
velocity  as  the  wave  propagates  along  the  symmetry  axis  towards  the  end  of  the  PDE 
configuration.  This  is  shown  in  Figure  47.  Note  that  the  shock  wave  passage  starts  out  at 
a  run  distance  of  23  mm  along  the  symmetry  axis  due  to  the  implosion  process  at  the 
mid-axis  inlet.  The  first  two  run  distance  points  are  also  discarded  due  to  the 
implementation  of  the  five  point  stencil  finite  difference  method  used  to  determine  the 
wave  velocity. 


Plot  of  Wave  Velocitv/CJ  Velocity  against  Run  Distance  along  Symmetry 

2  - 
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Figure  47,  Plot  of  wave  velocity/CJ  velocity  against  run  distance  along  symmetry 

axis 


The  snap  shot  pressure  and  temperature  contours  in  Figures  48  and  49  denoting 
the  point  just  after  impact  of  the  implosion  process  at  a  run  distance  of  23  mm  along  the 
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symmetry  axis  indicates  extremely  high  localized  pressure  and  temperature  regions  just 
aft  of  the  shock  front.  The  multi-dimensional  effects  of  the  wave  interactions  are 
explicitly  observed  in  the  temperature  contour.  The  pressure  and  temperature  history 
evaluated  along  the  symmetry  axis  in  Figures  50  and  51  allows  for  quantification  of  the 
pressure  and  temperature  states  arising  from  the  collision  of  the  two  implosion  shock 
waves.  The  localized  pressure  state  along  the  symmetry  axis  increased  to  as  high  as  6.4 
MPa  with  the  corresponding  temperature  at  3700  K.  In  comparison  to  the  characteristic 
CJ  pressure,  this  is  on  the  order  of  at  least  4  times  greater  and  the  detonation  wave  is 
highly  sustained  by  the  implosion. 
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Figure  48.  Pressure  contour  for  imploding  geometry  at  run  distance  of  23  mm 

along  symmetry  axis 


51 


E 


CO 


CO 

cc. 


0.1 


0.08 


0.06 


Symmetry  Length  (m) 


■ 

T 

3000 

2700 

1  800 

900 

Figure  49,  Temperature  contour  for  imploding  geometry  at  run  distance  of  23 

mm  along  symmetry  axis 


7,0E+06 

5.0E+06 

1 

i 

(Q 

i 

1 

a? 

i 

^  3.0E+06 

; 

1 

!  I 

CL 

1 

1 

. 

1 

1,0E+06 

' 

1 

O.OE+00  ^ 

0.02 

0.04 

0.06 

0.08 

Symmetry  Length  (m) 


Figure  50,  Pressure  history  of  imploding  wave  at  run  distance  of  23  mm  along 

symmetry  axis 
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Figure  51.  Temperature  history  of  imploding  wave  at  run  distance  of  23  mm 

along  symmetry  axis 

As  the  detonation  wave  front  propagates  along  the  symmetry  axis  to  a  run 
distanee  of  50  mm  as  shown  in  the  pressure  and  temperature  eontours  of  Figures  52  and 
53,  it  eatehes  up  with  the  spherieal  wave  front  extending  from  the  symmetry  length  to  the 
wall  of  the  tube  eonfiguration.  The  thin  region  behind  the  spherieal  wave  is  still  elearly 
being  driven  by  the  reaetion  proeess  as  it  stays  attaehed  to  the  shoek  front  above  C-J 
pressure  and  temperature  while  the  region  along  the  symmetry  axis  are  propagating  at 
extremely  high  overdriven  pressures  and  temperatures  fueled  by  the  eollision  of  the 
imploding  shoeks.  The  pressure  and  temperature  history  plots  along  the  symmetry  axis  in 
Figures  54  and  55  verifies  this  with  the  lead  shoek  front  at  a  pressure  and  temperature 
state  as  high  as  6.0  MPa  and  4000  K  respeetively. 
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Figure  52,  Pressure  contour  for  imploding  geometry  at  run  distance  of  50  mm 

along  symmetry  axis 
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Figure  53,  Temperature  contour  for  imploding  geometry  at  run  distance  of  50 

mm  along  symmetry  axis 
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Figure  55,  Temperature  history  of  imploding  wave  at  run  distance  of  50  mm 

along  symmetry  axis 


At  a  run  distance  of  98  mm  towards  the  end  of  the  tube,  the  pressure  and 
temperature  contours  of  Figures  56  and  57  depiets  the  merging  of  the  two  separate  waves 
into  an  almost  planar  wave  front  propagating  at  C-J  velocity.  The  localized  high  pressure 
region  stays  highly  attaehed  to  the  shock  front  followed  by  an  expansion  region.  A  triple 
point  discontinuity  can  be  observed  further  upstream  of  the  expansion  region  denoting  the 
interaction  of  the  longitudinal  and  transverse  wave  from  the  confined  walls  of  the  shoek 
tube.  The  pressure  and  temperature  states  along  the  symmetry  axis  pertaining  to  the 
leading  edge  of  the  shock  front  are  still  at  2.0  MPa  and  3000  K  respectively  from  Figures 
58  and  59  above  C-J  conditions,  which  highlight  the  fact  that  a  detonation  wave  still 
persists  and  will  eontinue  to  propagate  along  the  length  of  the  tube. 
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Figure  56,  Pressure  contour  for  imploding  geometry  at  run  distance  of  98  mm 

along  symmetry  axis 
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Figure  57,  Temperature  contour  for  imploding  geometry  at  run  distance  of  98 

mm  along  symmetry  axis 


Figure  58,  Pressure  history  of  imploding  wave  at  run  distance  of  98  mm  along 

symmetry  axis 
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Figure  59.  Temperature  history  of  imploding  wave  at  run  distance  of  98  mm 

along  symmetry  axis 


Analysis  of  the  implosion  process  along  both  the  mid-axis  of  the  inlet  and 
symmetry  axis  points  to  the  fact  that  a  detonation  shock  wave  persists  almost  throughout 
the  entire  geometry  and  is  sustained  both  by  the  imploding  wave  front.  This  highlights 
the  potential  of  introducing  a  toroidal  implosion  process  for  PDE  type  applications  and 
allows  for  larger  diameters  to  be  factored  into  such  configuration  designs  such  that 
shorter  tube  lengths  can  be  factored  in  without  the  need  to  account  for  DDT  run 
distances.  Subsequently,  a  parametric  study  will  be  conducted  to  evaluate  exactly  the 
extent  of  the  diameter  increase  such  that  a  detonation  shock  wave  can  be  sustained  for 
such  designs. 

D.  PARAMETRIC  STUDY 

Parametric  studies  were  performed  on  tube  diameters  of  114.3  mm  and  133.35 
mm  respectively  with  the  same  inlet,  reactant  and  boundary  conditions  specified  for  the 
toroidal  implosion  analysis.  Note  that  the  previous  tube  diameter  of  76.2  mm  from  the 
implosion  analysis  was  included  as  part  of  the  parametric  study.  Current  computational 
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resources  could  not  support  further  exploration  of  tube  diameters  beyond  133.35  mm  due 
to  huge  computational  demands  required  and  this  could  potentially  take  months  to 
compute  depending  on  the  test  case  model.  For  the  study,  geometry  of  the  shock  tube  as 
shown  in  Figure  18  remained  the  same  with  the  exception  of  the  diameter.  The  wave 
velocity  was  similarly  evaluated  along  the  mid-inlet  and  symmetry  axis  to  determine  if 
the  detonation  wave  can  be  sustained. 

1,  Parametric  Analysis  Along  Mid-inlet  Axis 

Figure  60  illustrates  the  detonation  wave  velocity  as  a  function  of  the  implosion 
shock  tube  radius  for  the  different  diameters  specified  along  the  mid- inlet  axis. 
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Figure  60,  Plot  of  wave  velocity/C J  Velocity  against  tube  radius  along  inlet  axis 
for  different  imploding  shock  tube  diameters 


As  the  tube  diameter  increases,  the  detonation  wave  velocity  starts  out  lower  as 

under-driven  detonations  due  to  decreased  confinement  effects.  An  analysis  of  the 

pressure  and  temperature  contours  for  diameters  of  114.3  mm  and  133.35  mm  as  the 
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imploding  wave  propagates  mid-way  through  the  radius  of  the  tube  shows  lower  pressure 
and  temperature  states  behind  the  lead  shoek  front  as  eompared  to  the  diameter  of  76.2 
mm.  Figures  61  and  62  illustrates  the  pressure  and  temperature  eontours  for  114.3  mm. 
Figures  63  and  63  illustrates  the  pressure  and  temperature  eontours  for  133.35  mm. 
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Figure  61.  Pressure  contour  for  imploding  geometry  at  run  distance  of  28.75  mm 
for  tube  diameter  of  114.30  mm  along  inlet  axis 
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Figure  62.  Temperature  contour  for  imploding  geometry  at  run  distance  of  28.75 
mm  for  tube  diameter  of  114.30  mm  along  inlet  axis 
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Figure  63.  Pressure  contour  for  imploding  geometry  at  run  distance  of  33.34  mm 
for  tube  diameter  of  133.35  mm  along  inlet  axis 
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Figure  64,  Temperature  contour  for  imploding  geometry  at  run  distance  of  33,34 
mm  for  tube  diameter  of  133,35  mm  along  inlet  axis 


For  the  tube  diameter  pertaining  to  114.30  mm,  the  phenomena  eorresponding  to 
the  coalescenee  of  tiny  spherical  wavelets  representing  localized  high  pressure  regions  at 
the  leading  edge  of  the  toroidal  wave  front  could  be  observed  similar  to  the  implosion 
analysis.  However  the  pressure  states  in  these  localized  regions  were  lower  in  comparison 
due  to  geometrical  differences  in  diameter.  Similarly  in  comparison  to  the  temperature 
contour  for  the  implosion  analysis,  there  is  a  thin  separation  region  between  the  reactant 
conditions  and  the  high  temperature  region  at  the  toroidal  front  with  an  increasing 
thickness  further  away  on  both  sides  of  the  mid-axis  inlet  indicating  a  decaying 
detonation  process  due  to  expansion  effects.  These  effects  are  however  greater  as 
compared  to  the  test  case  model  for  76.2  mm  used  for  the  implosion  analysis  and  this 
again  is  a  function  of  the  diameter  increase.  This  would  account  for  the  lower  detonation 
velocities  observed  as  the  tube  diameter  is  increased  from  76.2  mm  to  1 14.30  mm. 

For  the  case  study  of  133.35  mm,  the  pressure  contour  indicates  a  weak  decaying 
wave  propagating  mid-way  towards  the  symmetry  axis.  The  pressure  region  behind  the 
decaying  shock  front  is  almost  at  a  uniformly  low  pressure  with  a  magnitude  of 
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approximately  half  of  CJ  pressure.  However  due  to  the  geometry  of  the  imploding  tube,  a 
seeond  refleetion  wave  from  the  head  end  of  the  tube  ean  be  observed  interaeting  in  a 
direetion  perpendieular  to  the  weak  wave  propagation  along  the  mid-inlet  axis.  This 
refleetion  wave  represents  a  strong  detonation  wave  with  loealized  pressure  regions  just 
behind  the  toroidal  wave  front  starting  from  the  head  end  exeeeding  two  times  of  C-J 
pressure.  Likewise,  as  for  the  temperature  eontour,  there  exists  a  thiek  uniform  region  of 
separation  between  the  reaetant  temperature  and  high  temperature  region  along  the 
toroidal  wave  front  with  an  interaeting  refleetion  wave  from  the  shoek  tube  head  end.  The 
weaker  pressure  and  temperature  states  behind  the  deeaying  shoek  front  aeeount  for  even 
lower  detonation  veloeities  as  eompared  to  the  tube  diameter  eorresponding  to  114.30 
mm. 

As  the  wave  implodes  towards  the  symmetry  axis,  the  proeess  fueling  the  lead 
shoek  front  is  not  supported  and  eontinues  to  weaken  and  deeay  due  to  expansion  effeets 
with  the  exeeption  of  the  test  ease  pertaining  to  133.35  mm.  The  seeond  refleeted  wave 
from  the  head  end  along  the  symmetry  axis  reaehes  interaets  with  the  transverse  wave 
propagation  along  the  mid-inlet  axis,  whieh  reinforees  the  lead  shoek  front  and  results  in 
a  sharp  inerease  in  detonation  veloeity  towards  0.8  of  C-J  veloeity  as  shown  in  Figure  60. 
The  respeetive  pressure  and  temperature  eontours  in  Figures  65  to  66  and  67  to  68 
eorresponding  to  the  test  ease  of  1 14.30  mm  and  133.35  mm  illustrate  the  observed  wave 
propagation  behavior  as  the  shoek  front  implodes  towards  the  symmetry  axis. 
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Figure  65.  Pressure  contour  for  imploding  geometry  at  run  distance  of  52  mm 
for  tube  diameter  of  114.30  mm  along  inlet  axis 
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Figure  66.  Temperature  contour  for  imploding  geometry  at  run  distance  of  52 
mm  for  tube  diameter  of  114.30  mm  along  inlet  axis 
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Figure  67.  Pressure  contour  for  imploding  geometry  at  run  distance  of  62  mm 
for  tube  diameter  of  133.35  mm  along  inlet  axis 
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Figure  68.  Temperature  contour  for  imploding  geometry  at  run  distance  of  62 
mm  for  tube  diameter  of  133.35  mm  along  inlet  axis 
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2,  Parametric  Analysis  Along  Symmetry  Axis 

Along  the  symmetry  axis,  the  detonation  proeess  is  supported  throughout  as  ail  3 
test  eases  pertaining  to  different  tube  diameter  geometries  is  observed  to  start  from  an 
overdriven  detonation  wave  before  decaying  to  a  steady  state  wave  propagating  towards 
C-J  velocity.  A  plot  of  the  normalized  wave  velocity  to  C-J  velocity  against  run  distance 
along  symmetry  axis  in  Figure  69  depicts  the  observed  wave  propagation  behavior. 


Plot  of  Wave  Velocitv/CJ  Velocity  against  Run  Distance  along  Symmetry 

Axis 
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Figure  69,  Plot  of  wave  velocity/CJ  Velocity  against  run  distance  along  symmetry 
axis  for  different  imploding  shock  tube  diameters 


The  test  case  pertaining  to  tube  diameter  of  1 14.3  mm  starts  out  immediately  at  60 
mm  along  the  symmetry  axis  as  the  geometry  of  the  wave  front  along  the  radial  direction 
is  such  that  at  the  point  of  implosion  is  almost  planar  throughout  till  the  start  point.  The 
implosion  process  along  the  wave  propagation  drives  the  pressure  and  temperature  along 
the  symmetry  plane  to  exceedingly  high  values  to  almost  50  times  C-J  pressure  as  shown 
in  the  pressure  and  temperature  history  of  the  starting  wave  at  run  distance  of  60  mm  as 
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illustrated  in  Figures  70  and  71.  This  eauses  a  sharp  inerease  in  C-J  veloeity  to  almost  1.7 
of  C-J  veloeity,  whieh  eauses  the  detonation  wave  to  start  out  overdriven.  There  is  also  a 
second  reflected  wave  from  the  head  end  of  the  tube,  which  is  propagating  behind  the 
starting  wave  as  shown  in  the  pressure  contour  in  Figure  72. 
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Figure  70,  Pressure  history  of  imploding  wave  at  run  distance  of  60  mm  along 

symmetry  axis  (diameter  114,3  mm) 
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Figure  71,  Temperature  history  of  imploding  wave  at  run  distance  of  60  mm 
along  symmetry  axis  (diameter  114,3  mm) 
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Figure  72.  Pressure  contour  for  imploding  geometry  at  run  distance  of  60  mm 
for  tube  diameter  of  114.3  mm  along  symmetry  axis 
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The  detonation  proeess  is  supported  throughout  the  symmetry  axis  as  both 
pressure  and  temperature  states  deeay  towards  C-J  eonditions.  Even  as  the  wave  approaeh 
towards  the  end  of  the  tube  at  a  run  distance  of  91  mm,  the  pressure  and  temperature 
states  at  the  1-D  lead  shock  front  along  the  symmetry  axis  as  denoted  by  the  respective 
temperature  and  pressure  histories  in  Figures  73  and  74  are  still  above  C-J  conditions. 
This  clearly  indicates  that  the  detonation  propagation  process  is  still  supported.  In  2-D, 
the  pressure  contour  in  Figure  75  at  a  run  distance  of  91  mm  illustrates  a  non-planar 
shock  front  followed  by  a  closely  coupled  high  pressure  region.  In  this  region  along  the 
transverse  direction,  the  pressure  decreases  further  away  from  the  symmetry  axis  before 
increasing  again  towards  the  wall  of  the  shock  tube.  Although  there  is  a  decrease  of 
pressure  in  the  transverse  direction,  the  states  in  this  closely  coupled  are  still  above  C-J 
conditions,  which  highlights  the  potential  of  increasing  the  diameter  of  the  tube  even 
further  and  still  maintain  sustenance  of  the  detonation  wave. 


Figure  73.  Pressure  history  of  imploding  wave  at  run  distance  of  91  mm  along 

symmetry  axis  (diameter  114,3  mm) 
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Figure  74,  Temperature  history  of  imploding  wave  at  run  distance  of  91  mm 
along  symmetry  axis  (diameter  114,3  mm) 
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Figure  75.  Pressure  contour  for  imploding  geometry  at  run  distance  of  91  mm 
for  tube  diameter  of  114.3  mm  along  symmetry  axis 
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For  the  test  ease  pertaining  to  the  tube  diameter  of  133.35  mm,  there  are  2 
meehanisms  sustaining  the  detonation  wave  propagation.  The  first  meehanism  is  due  to 
the  refleetion  of  the  seeondary  wave  from  the  head  end  of  the  shoek  tube  while  the 
seeond  meehanism  is  due  to  the  implosion  proeess  itself.  At  implosion,  the  wave  starts 
out  overdriven  with  pressure  and  temperature  states  behind  the  lead  shoek  front  way 
above  C-J  eonditions.  For  the  speeified  geometry,  the  seeondary  refieeted  wave  merges 
with  the  transverse  imploded  wave,  whieh  magnifies  the  pressures  and  temperatures 
behind  the  lead  shoek  front.  This  ean  be  observed  in  the  pressure  and  temperature 
histories  of  the  wave  propagation  evaluated  along  the  symmetry  axis  at  a  run  distanee  of 
30  mm  in  Figures  76  and  77,  whieh  reaehes  32.5  MPa  (~20  times  C-J  pressure)  and  5000 
K  (~  1.7  times  C-J  temperature).  The  same  2-D  pressure  contour  in  Figure  78  illustrates 
an  observed  two  distinct  wave  structure  phenomena,  a  detonation  wave  that  is 
propagating  at  super  C-J  velocity  and  a  weak  decaying  wave  that  is  remnant  of  the 
previous  toroidal  wave  front  that  is  slowing  down  due  to  expansion  effects. 


Figure  76,  Pressure  history  of  imploding  wave  at  run  distance  of  30  mm  along 

symmetry  axis  (diameter  133,35  mm) 
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Figure  77,  Temperature  history  of  imploding  wave  at  run  distance  of  30  mm 
along  symmetry  axis  (diameter  133,35  mm) 
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Figure  78.  Pressure  contour  for  imploding  geometry  at  run  distance  of  30  mm 
for  tube  diameter  of  133.35  mm  along  symmetry  axis 
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The  super  detonation  wave  travelling  along  the  symmetry  axis  will  eventually 
catch  up  and  overtake  the  weak  expansion  toroidal  wave  front.  This  is  illustrated  in  the 
pressure  contour  in  Figure  79  as  the  wave  approaches  the  end  of  the  shock  tube  length  at 
a  mn  distance  of  88  mm.  The  super  detonation  wave  also  expands  towards  the  walls 
eventually  replacing  the  weak  expansion  wave  front  as  it  propagates  further  along  the 
symmetry  axis. 
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Figure  79,  Pressure  contour  for  imploding  geometry  at  run  distance  of  88  mm 
for  tube  diameter  of  133,35  mm  along  symmetry  axis 


The  1-D  pressure  and  temperature  histories  of  the  same  wave  along  the  symmetry 
axis  in  Figures  80  and  81  still  shows  the  detonation  wave  being  supported  even  as  it 
propagates  towards  the  tube  end.  Here  the  C-J  pressure  and  temperature  states  just  behind 
the  lead  shock  front  are  still  at  4.0  MPa(~  2.5  C-J  pressure)  and  3200  K  (~  1.1  C-J 
temperature),  which  will  eventually  decay  to  steady-state  CJ  velocity. 
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Figure  81.  Temperature  history  of  imploding  wave  at  run  distance  of  88 
along  symmetry  axis  (diameter  133,35  mm) 


The  parametric  study  highlights  the  immense  potential  of  utilizing  a  shock 
focusing  technique  such  as  toroidal  implosion  to  increase  specific  thrust  in  pulse 
detonation  engines  with  no  requirements  to  factor  in  tube  length  for  DDT  transitions  as 
the  diameter  of  the  tube  is  further  increased.  Even  when  the  diameter  was  increased  to 
beyond  1.5  times  the  original  length,  the  detonation  process  was  still  supported 
throughout  as  shown  by  the  simulation  results.  The  mechanisms  supporting  the 
detonation  process  was  due  to  the  implosion  process  itself  as  well  as  the  secondary 
reflection  wave  from  the  tube  head  end.  It  would  have  been  ideal  to  obtain  the  optimal 
tube  diameter  at  which  the  implosion  process  fails  to  support  the  detonation  process; 
however  it  was  not  possible  to  go  beyond  the  tube  diameter  of  133.35  mm  due  to  the 
limitation  of  the  available  computational  resources.  A  two-fold  increase  in  tube  diameter 
results  in  a  four-fold  increase  in  required  computational  resource  due  to  the  number  of 
mesh  size.  It  is  evident  from  the  outcome  of  the  parametric  study  that  diameters  way 
beyond  133.35  mm  is  possible  to  support  the  detonation  process  due  to  the  simultaneous 
interaction  and  reinforcement  of  not  only  the  implosion  process  itself,  but  also 
supplemented  by  the  secondary  reflection  wave  from  the  head  end  of  the  tube. 
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V.  WAY  FORWARD  AND  RECOMMENDATIONS 


Steps  to  conduct  further  parametric  studies  to  determine  the  optimal  diameter 
such  that  detonation  process  can  be  supported  will  be  addressed  in  future  endeavors.  The 
specialized  CFD  solver  from  Metacomp  Technologies  requires  extensive  parallel  node- 
to-node  communications  based  on  the  computational  cell  size  required  to  capture  the 
complex  physics  governing  the  shock  wave  propagation  process.  This  results  in  latency 
issues  that  contribute  to  substantial  computational  downtime  whereby  current 
supercomputing  architecture  on  campus  using  infiniband  interconnects  is  unable  to 
support.  Subsequent  studies  will  be  ported  over  to  a  high  performance  computing 
architecture  off-campus  with  minimal  latency  and  as  such  simulation  time  can  be 
drastically  reduced  in  order  to  allow  for  reasonable  time  for  analysis  of  results. 

Additionally,  it  is  desired  not  only  to  maximize  the  tube  diameter  based  on  the  given 
geometry,  but  also  the  location  of  the  inlet  for  the  starting  toroidal  wave  front.  The  location 
of  the  inlet  from  the  head  tube  end  will  not  only  optimize  the  imploding  detonation  wave 
towards  the  symmetry  axis  but  also  optimizes  the  secondary  reflected  wave  as  it  catches  up 
with  the  toroidal  wave  front  and  aids  as  a  supplementary  mechanism  to  support  the 
detonation  wave  propagation  throughout  the  shock  tube  along  the  symmetry  axis. 

Once  the  optimal  configuration  geometry  has  been  established  for  imploding 
detonations  for  hydrogen-air,  future  endeavors  will  also  address  test  cases  for 
hydrocarbon  fuels  such  as  ethylene-air  and  JP  10-air.  Moreover,  optimal  tube  geometries 
can  be  further  determined  by  using  the  practical  computational  cell  size  established  for 
this  thesis  study  as  well  as  a  set  of  viable  reduced  set  of  chemical  kinetic  reaction  models 
specific  to  these  fuel-oxidizer  mixtures.  Note  that  JPIO  is  a  military  qualified  jet  aviation 
fuel  that  due  to  its  high  energy  density  represents  the  best  realizable  fuel  of  choice  for 
PDE  applications. 

The  future  should  focus  on  building  an  experimental  test  rig  based  on  the 
imploding  geometry  to  not  only  validate  the  simulation  results,  but  in  the  process  allow 
for  a  paradigm  shift  in  future  PDE  designs  that  need  not  factor  in  multiple  tube 
configurations  and  as  such  maximizing  specific  thrust. 
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VI.  CONCLUSION 


A  novel  method  to  increase  specific  thrust  in  PDEs  via  imploding  detonations  has 
been  investigated  in  this  study.  High  fidelity  simulations  were  required  in  order  to  capture 
the  physics  governing  such  complex  flow  fields.  In  this  study,  a  practical  computational 
mesh  size  on  the  order  of  one-quarter  induction  length  was  determined  and  deemed 
sufficient  based  on  a  reduced  set  of  18  reactions  and  9  species  chemical  reaction  model 
for  hydrogen-air  mixture.  It  was  also  found  that  the  corresponding  computed  detonation 
cell  width  coincided  with  actual  experimental  data  for  the  same  fuel-oxidizer  mixture. 

A  set  of  simulations  were  carried  out  on  an  enhanced  PDE  design  configuration 
utilizing  implosion  concepts  and  results  indicate  that  the  detonation  wave  was  supported 
with  pressure  and  temperature  states  persisting  above  CJ  conditions  behind  the  lead  shock 
front  throughout  the  entire  geometry,  further  parametric  studies  were  conducted  to 
evaluate  and  determine  the  extent  of  diameter  increase  for  such  design  configurations  and 
it  was  found  that  the  detonation  propagation  process  was  still  well  supported  even  up  to 
1.75  times  the  original  diameter.  This  represents  a  three-fold  increase  in  thrust  per  PDE 
tube,  which  roughly  translates  into  a  two-fold  increase  in  specific  thrust  for  similar  total 
thrust  PDE  systems.  However,  there  is  a  need  to  evaluate  much  bigger  diameters.  Due  to 
the  limitations  of  current  parallel  computing  architecture,  future  endeavors  will  be 
addressed  by  using  the  Department  of  Defense  (DoD)  High  Performance  Computing 
systems  that  consists  of  a  parallel  computing  platform  with  minimal  latency  in  order 
determine  the  maximum  diameter  at  which  the  detonation  structure  fails  to  persist  in  the 
specified  geometry. 
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APPENDIX  A:  NASA  CHEMICAL  EQUILIBRIUM  ANALYSIS 

OUTPUT 


NASA-GLENN  CHEMICAL  EQUILIBRIUM  PROGRAM  CEA2 ,  FEBRUARY  5,  2004 
BY  BONNIE  MCBRIDE  AND  SANFORD  GORDON 
REFS:  NASA  RP-1311,  PART  I,  1994  AND  NASA  RP-1311,  PART  II,  1996 

■k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k'k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k 

'k'k'k'k'k'k'k'k 


output  massf 
prob  case=00002142  det 
phi=l 

p (atm) =1 
output  short 
output  trace=le-10 

reac 

oxid  Air  t,k=  298  wt%=  100. 
fuel  H2  t,k=  298  wt%=  100. 
end 

DETONATION  PROPERTIES  OF  AN  IDEAL  REACTING  GAS 

CASE  =  00002142 

REACTANT 

TEMP 
K 

OXIDANT  Air 

298 . 000 

FUEL  H2 

298.000 

0/F=  34.29623  %FUEL=  2.833164  R,EQ.RATIO=  1.000000  PHI , EQ .  RATIO= 

1 . 000000 

UNBURNED  GAS 


PI, 

BAR 

1 . 0132 

Tl, 

K 

298.00 

HI, 

KJ/KG 

-4.36 

Ml, 

(1/n) 

21.008 

GAMMA 1 

1.4015 

SON 

VELl, M/SEC 

406.6 

BURNED  GAS 


WT  FRACTION 
(SEE  NOTE) 

1 . 0000000 
1 . 0000000 


ENERGY 

KJ/KG-MOL 

-129.895 

-4.325 
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P,  BAR 

15.801 

T,  K 

2943.63 

RHO,  KG/CU  M 

1.5500  0 

H,  KJ/KG 

1333.33 

U,  KJ/KG 

313.90 

G,  KJ/KG 

-29833.7 

S,  KJ/ (KG) (K) 

10.5880 

M,  (1/n) 

24.008 

(dLV/dLP) t 

-1 . 00954 

(dLV/dLT) p 

1.2062 

Cp,  KJ/ (KG) (K) 

3.3542 

GAMMAS 

1 .1637 

SON  VEL, M/SEC 

1089.2 

DETONATION  PARAMETERS 

P/Pl 

15.595 

T/Tl 

9.878 

M/Ml 

1 . 1428 

RHO/ RHO 1 

1 . 8042 

DET  MACH  NUMBER 

4 . 8335 

DET  VEL, M/SEC 

1965.1 

MASS  FRACTIONS 

*Ar 

1.2550-2 

*CO 

1.3093-4 

*C02 

2 . 6525-4 

COOH 

2 . 6149-9 

*H 

2.4760-4 

HCO 

1.0225-9 

HNCO 

2.395-10 

HNO 

3.1235-6 

HN02 

8 . 9839-7 

HN03 

1.028-10 

H02 

1.5250-5 

*H2 

2 . 6254-3 

HCOOH 

3.735-10 

H20 

2.2035-1 

H202 

2.3080-6 

*N 

1 .1146-6 

*NH 

4.5040-7 

NH2 

2.7000-7 

NH3 

3.4969-7 

NH20H 

4.217-10 

*NO 

9.4235-3 

N02 

6.6339-6 

*N2 

7.2938-1 

N20 

2.7711-6 

*0 

1.3571-3 

*0H 

1.3522-2 

*02 

1.0117-2 

03 

6.2078-9 
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*  THERMODYNAMIC  PROPERTIES  FITTED  TO  20000. K 


NOTE.  WEIGHT  FRACTION  OF  FUEL  IN  TOTAL  FUELS  AND  OF  OXIDANT  IN  TOTAL 
OXIDANTS 
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APPENDIX  B:  CONVERGENCE  PLOTS  FOR  RESIDUAL  ERRORS 


A,  GLOBAL  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.25  MM 


B.  INNER  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.25  MM 
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c 


GLOBAL  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.50  MM 


k  Applications  Places  System 


)e+04  4.22S1e-07 


■  ▼elixee  of  tfie  RHS  ooIxn 


MAXIMDK  KENniniM 


S . B9&9e*03  5 . 9431e< 


All  residuals  whose  maximum  ks  below:  jl  .0e-l4  will  not  be  plotted  Apply  | 
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D.  INNER  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.50  MM 
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E, 


GLOBAL  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  1.00  MM 


F.  INNER  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  1.00  MM 
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G. 


GLOBAL  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.0625  MM 
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I 


GLOBAL  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.05  MM 


J.  INNER  RESIDUAL  ERROR  PLOT  FOR  MESH  SIZE  OF  0.05  MM 
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K 


GLOBAL  RESIDUAL  ERROR  PLOT  FOR  IMPLOSION  (76.2  MM) 


L.  INNER  RESIDUAL  ERROR  PLOT  FOR  IMPLOSION  (76.2  MM) 
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M.  GLOBAL  RESIDUAL  ERROR  PLOT  FOR  IMPLOSION  (1 14.30  MM) 


N.  INNER  RESIDUAL  ERROR  PLOT  FOR  IMPLOSION  (1 14.30  MM) 
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o 


GLOBAL  RESIDUAL  ERROR  PLOT  FOR  IMPLOSION  (133.35  MM) 


P.  INNER  RESIDUAL  ERROR  PLOT  FOR  IMPLOSION  (133.35  MM) 
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